This workshop will provide you with an introduction to manipulating raster spatial data using the R package terra. terra and its predecessor, raster are widely used for spatial data manipulation in R. No prior experience of spatial data is assumed, but this short workshop will not have time to delve into some important aspects of spatial data such as projections.

terra’s predecessor raster had many of the same functions as terra, and I will mention how functions have changed or been renamed which might be helpful for people migrating from using raster. In the words of the creator of both packages: terra is simpler, faster and can do more, so definitely switch to terra if you are still using raster!

raster users note: Boxes like this will highlight differences between the raster package and terra.

Resources:

The official terra tutorial page is https://rspatial.org/spatial/index.html
An excellent reference to all the functions in terra is https://rspatial.github.io/terra/reference/terra-package.html
This website is an good guide to doing GIS work in R using both terra and sf packages: https://mgimond.github.io/Spatial/index.html

Prerequisites

You will need the terra package installed, which can be done by running install.packages("terra"). If you have problems, there are more details about installing it here. We can then load it

library(terra)
#> terra 1.7.71

Spatial data

There are basically two types of spatial data: vector and raster

Vector data

Can be points, lines or polygons. Useful for representing things like survey locations, rivers, and boundaries.

pts <- rbind(c(3.2,4), c(3,4.6), c(3.8,4.4), c(3.5,3.8), c(3.4,3.6), c(3.9,4.5)) |>
  vect()

lnes <- as.lines(vect(rbind(c(3,4.6), c(3.2,4), c(3.5,3.8)))) |>
  rbind(as.lines(vect(rbind(c(3.9, 4.5), c(3.8, 4.4), c(3.5,3.8), c(3.4,3.6)))))

lux <- vect(system.file("ex/lux.shp", package = "terra"))

par(mfrow = c(1,3))

plot(pts, axes = F, main = "Points")
plot(lnes, col = "blue", axes = F, main = "Lines")
plot(lux, "NAME_2", col = terrain.colors(12), las = 1, axes = F, main = "Polygons")


par(mfrow = c(1,3))

Raster data

Raster data is a grid of rectangles, normally called cells. Each cell has a value, making rasters useful for storing continuous data, such as temperature and elevation.

Here is an example of raster data, where each cell in the raster represents elevation.

elev <- system.file("ex/elev.tif", package = "terra") |>
  rast() |>
  aggregate(fact = 2)

plot(elev, las = 1, main = "Elevation map")

elev |>
  as.polygons(aggregate = FALSE, na.rm = FALSE) |>
  lines(col = "grey40", lwd = 0.2)

Getting started

Lets start by creating our own raster. We can create rasters from scratch and load them from a file using the function rast(). We can create a simple raster by specifying the x and y limits for the raster and the resolution (how big each cell is).

raster users note: rast() replaces raster()

#create raster
ras <- rast(xmin = 0, xmax = 10, ymin = 0, ymax = 10, resolution = 2)

#see what we've created
ras
#> class       : SpatRaster 
#> dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
#> resolution  : 2, 2  (x, y)
#> extent      : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84

The figure below shows what most of the terms above refer to. As you can see, you don’t need to use all the terms to define a raster. Couple of other points:

Every object in R has a class, such as data.frame and as you can see, rasters in terra are of class SpatRaster.
We did not tell rast() which coordinate reference system to use, so it defaults to using longitude latitude coordinates, also known as EPSG 4326.

raster users note: rasters are now SpatRaster class not RasterLayer

But what does the raster we created actually look like when plotted. Lets see. All we need is plot()

plot(ras)

Why is there no plot? Because the raster we created is empty; there are no values associated with the the cells. Lets assign some values to each cell in the raster and try again. First we will find out how many cells are in our raster using `ncell()

ncell(ras)
#> [1] 25

Ok, now we know this lets give our raster cells values from 1 to 25:

values(ras) <- 1:25

plot(ras)

Now our raster has values, we get a plot! Each cell has an integer value between 1 and 25, with cell values increasing from left to right and top to bottom. So the values start being “filled up” in the top left, and finish in the bottom right.

Lets have another look at our raster properties

ras
#> class       : SpatRaster 
#> dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
#> resolution  : 2, 2  (x, y)
#> extent      : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source(s)   : memory
#> name        : lyr.1 
#> min value   :     1 
#> max value   :    25

We can now see a few extra pieces of information compared to last time:

sources(s): where is the data held on your computer? It says memory for this raster, indicating that the raster is in the computer memory. Rasters can also be held on your hard disk, in which case this will be the file name of the raster. We won’t go into details here, but terra is smart about loading data into memory, only doing so when it needs to and it thinks it will have enough space.
name: what is the raster called?
min value & max value: the minimum and maximum values in the raster

Ok, now we understand the basic structure of a raster, lets load some data.

Exploring some real world data

Sea surface temperature can be measured by satellites, and the National Oceanic and Atmospheric Administration (NOAA) in the U.S. has been doing so since 1972! Daily sea surface temperature for the world from 1985 to the present is available via the NOAA Coral Reef Watch website.

We are going to explore sea surface temperature (SST) data for the Great Barrier Reef in Australia. All the data is in the data folder in the Github repository, and full details on how I got the data are in the data_prep.R script in that folder.

Load the data

We can use the same rast() command that we used to make our own raster to load data The data is stored in GeoTiff (.tif) format, which is widely used for raster data. You can get a full list of data formats that terra can read and write by running gdal(drivers = TRUE).

sst <- rast("../data/gbr_temp_2023_05_31.tif") #load the data

sst #view the properties
#> class       : SpatRaster 
#> dimensions  : 801, 1101, 1  (nrow, ncol, nlyr)
#> resolution  : 0.04999999, 0.05  (x, y)
#> extent      : 110, 165.05, -45, -4.95  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : gbr_temp_2023_05_31.tif 
#> name        : CRW_SST_2 
#> min value   :      9.59 
#> max value   :     30.46 
#> unit        :   Celsius

If the data loaded correctly, you should see the properties as shown above. We have all the properties I described earlier, but there are a couple of things worth noting:

the coordinate reference system (crs) is in lon/ lat, which means the unit of measurement is degrees. So the extent coordinates are in degrees, i.e. xmin = 110\(^\circ\), xmax = 165.05\(^\circ\), ymin = -45.00\(^\circ\) and ymax = 165.05\(^\circ\), and resolution is also in degrees, i.e 0.05\(^\circ\).
the source is the filename of the data we loaded. This means it is on disk, not in memory. We can double check by running inMemory(sst) which should return FALSE.
unit is Celsius. Most rasters don’t come with the units, so you will normally need to look at metadata to find this information.

Plot the data

We want to plot our data to see what it looks like. This is the first thing you should do with almost any spatial data to do a quick check that it looks right

plot(sst)

The values in the legend seem about what we would expect for water temperatures, ranging from ~10 - 30\(^\circ\) Celsius.

The white areas of the plot are cells with NA values, which in this case are land; there are no sea surface temperature data for the land! We can see the land mass of Australia taking up much of the map and Papua New Guinea to the north.

We can set the NA values to a different colour if we want, which can be helpful if you want to see some NA cells that are getting lost against the other colours.

plot(sst, colNA = "black")

Note that every cell within our raster has to have a value.

Load and plot some vector data

To put our raster data in context, lets plot the Great Barrier Reef Marine Park boundary (dowloaded from here). This is vector data. You might be used to handling vector data with the sf package, but terra can also be used for vector data manipulation. We can load a vector using vect()

raster users note: terra has its own methods for handling vector data unlike raster which used the sp package for vector data handling. Vector data in terra are SpatVector objects; different from sf objects.

#load the park boundary
gbr_boundary <- vect("../data/gbr_MPA_boundary.gpkg")

We can plot vector boundaries on top of raster data using lines():

#plot the SST data and the boundary on top
plot(sst)
lines(gbr_boundary) #this plots the vector lines over the raster data

We now see the outline of the Great Barrier Reef marine park on top of the sea surface temperature data we plotted before.

Crop and mask data

The data we have at the moment is for a much larger area than just the Great Barrier Reef marine park. To get only that data we need to crop and mask the raster data using the marine park boundary.

Cropping means that we keep only the data inside the extent of the vector we are using. Mask means that all the data outside the vector is set to NA or some other value we specify. Lets have a look how this works.

First lets have a look at the extent of the marine park boundary. We can get the extent of a raster or vector using ext(). We need to convert this into a SpatVector object for plotting using vect(). We only need to do this for plotting; when we crop, we can just use the marine park boundary as the input.

raster users note: ext() replaces extent()

gbr_boundary_extent <- ext(gbr_boundary) |>
  vect()

plot(sst)
lines(gbr_boundary)
lines(gbr_boundary_extent, col = "blue")

Cropping means we remove everything outside the extent (blue box) of our polygon. Masking sets all values outside our polygon to NA.

So when we crop, we get only the area within the blue box.

We crop using the crop() function, using the raster we want to crop as the first argument and the vector we are cropping with second.

#crop
sst_cropped <- crop(sst, gbr_boundary)

#plot
plot(sst_cropped)
lines(gbr_boundary)

Now we have cropped our raster, we can mask it so that we only have values for the area within the marine park boundary. We do this using mask:

#mask
sst_cropped_masked <- mask(sst_cropped, gbr_boundary)

#plot
plot(sst_cropped_masked)
lines(gbr_boundary)

Now we only see raster values for cells that are within the marine park boundary. But remember that the areas that are white, still have values, they are just NA values.

Often we want to crop and mask one after the other, and you can do this in one command using crop(sst, gbr_boundary, mask = TRUE).

For reference, here is a figure comparing what crop, mask and crop(mask = TRUE) do:

par(mfrow = c(2,2))

plot(sst, main = "Original raster")
lines(gbr_boundary)

plot(sst_cropped, main = "Cropped")
lines(gbr_boundary)

sst |>
  mask(gbr_boundary) |>
  plot(main = "Masked")
lines(gbr_boundary)

plot(crop(sst, gbr_boundary, mask = TRUE), main = "Cropped and masked")
lines(gbr_boundary)


par(mfrow = c(1,1))

Why not just mask rather than crop and mask? As we see in the figure above, this would mean we have a lot of area we are not interested in and even though most of those cells would be NA they take up space in our raster, so it is not efficient.

Raster values

Now we have cropped and masked our original raster to get only data within the area we are interested in, we can start exploring the values. terra has several functions that can help us do this easily.

Histogram

We can get a histogram of all the values in our raster using the hist() function, which is equivalent to the base R function.

hist(sst_cropped_masked)

The x-axis is in the units of our raster values; temperature in degrees Celsius in our case. The y-axis is frequency; how many cells in our raster have those values.

You can change the histogram using the same arguments you use with hist() in base R. For example, lets increase the number of bars we have:

hist(sst_cropped_masked, breaks = 100)

Frequency table

We can get a frequency table of values in our raster using freq(). This is essentially the same information that is shown in the histogram in graphical format.

freq(sst_cropped_masked)
#>   layer value count
#> 1     1    20    31
#> 2     1    21   223
#> 3     1    22   396
#> 4     1    23   740
#> 5     1    24  2861
#> 6     1    25  2827
#> 7     1    26  2194
#> 8     1    27  2638
#> 9     1    28   444

The default is to round values to the nearest integer. To get more integers we can do:

freq(sst_cropped_masked, digits = 1) |>
  head() #this is a long table: just show the first few values
#>   layer value count
#> 1     1  19.8     1
#> 2     1  19.9     3
#> 3     1  20.0     6
#> 4     1  20.1     3
#> 5     1  20.2     3
#> 6     1  20.3     7

We can also check how many NA values are in our raster:

freq(sst_cropped_masked, value = NA)
#>   layer value count
#> 1     1    NA 50850

Statistics

We can use the same summary() command that is used in base R to get a summary of the statistical information for an entire raster.

summary(sst_cropped_masked)
#>    CRW_SST_2    
#>  Min.   :19.82  
#>  1st Qu.:24.25  
#>  Median :25.05  
#>  Mean   :25.19  
#>  3rd Qu.:26.49  
#>  Max.   :27.74  
#>  NA's   :50850

These values are, for example, the mean value of all raster values excluding NAs.

To get individual statistical values, we need to use global(). For example to get the mean value of a raster:

global(sst_cropped_masked, "mean")
#>           mean
#> CRW_SST_2  NaN

Hmmm, what went wrong? The default is for global() to include all values in the raster, and since we have lots of NAs, the result is NaN. We need to set na.rm = TRUE to exclude the NA values.

global(sst_cropped_masked, "mean", na.rm = TRUE)
#>               mean
#> CRW_SST_2 25.18745

The object returned by global() is a data frame, so if you want just the value, you need to do:

sst_mean <- as.numeric(global(sst_cropped_masked, "mean", na.rm = TRUE))

sst_mean
#> [1] 25.18745

raster users note: global() replaces cellStats(), and the default is na.rm = FALSE, whereas the default for cellStats() was na.rm = TRUE.

Classifying

We might want to break our raster values into groups. For example, we could say that all temperatures that are below the mean temperature are classified as “cooler” and all temperatures above the mean are “warmer”. We can do this using the classify function.

raster users note: classify() replaces reclassify()

#first we create a matrix that will be used for the classification
# all values >= 0 and less than the mean become 1
# all values greater than the mean become 2
reclass_matrix <- c(0, sst_mean, 1,
                    sst_mean, Inf, 2) |>
  matrix(ncol = 3, byrow = TRUE)

#now we classify our raster using this matrix
sst_reclassed <- classify(sst_cropped_masked, reclass_matrix)

#plot the result
plot(sst_reclassed)

A gut check tells us this looks right; the warmer areas are in the north, nearer to the equator.

This plot is ok, but it would be better if the colours were more appropriate and the legend gave some useful information.

plot(sst_reclassed, col = c("blue", "red"), plg = list(legend = c("Cooler", "Warmer")), las = 1) #the las = 1 argument just rotates the y-axis labels so that they are horizontal

Raster math

The great thing about rasters are you can do maths with them! For example, doing sst_reclassed + 1 just adds one to each raster value, and doing sst_reclassed*2 multiplies each raster value by two.

As an example, lets convert our temperature data into Fahrenheit for our confused colleagues in the U.S. The conversion from Celsius to Fahrenheit is: Fahrenheit = (Celsius * 1.8) + 32.

#do the conversion
sst_fahrenheit <- (sst_cropped_masked*1.8) + 32

#plot our new raster
plot(sst_fahrenheit)

Nice maps

Lets also make the map a bit prettier, choosing more appropriate colours and adding a scale bar. There are huge number of options for plotting, see the help file ?help for details. There are many great plotting packages such as tmap and ggplot that can be used to make maps, but I will not cover those here. The Geocomputation with R website is an excellent resource on map making and geospatial data in R in general.

plot(sst_fahrenheit, col = hcl.colors(50, palette = "RdYlBu", rev = TRUE), las = 1) #see ?hcl.colors for more info on colour palettes. rev= TRUE because we want to reverse the palette: red colours are the highest values and blue the lowest
lines(gbr_boundary)
sbar(d = 400, type = "bar", divs = 2, below = "km") #400km scale bar with 2 divisions and "km" written below

Click to see a range of colour palettes

The swatches below show all the colour scales available via hcl.colors(). The code is directly from ?hcl.colors().

## color swatches for HCL palettes
hcl.swatch <- function(type = NULL, n = 5, nrow = 11,
  border = if (n < 15) "black" else NA) {
    palette <- hcl.pals(type)
    cols <- sapply(palette, hcl.colors, n = n)
    ncol <- ncol(cols)
    nswatch <- min(ncol, nrow)

    par(mar = rep(0.1, 4),
        mfrow = c(1, min(5, ceiling(ncol/nrow))),
        pin = c(1, 0.5 * nswatch),
        cex = 0.7)

    while (length(palette)) {
        subset <- 1:min(nrow, ncol(cols))
        plot.new()
        plot.window(c(0, n), c(0, nrow + 1))
        text(0, rev(subset) + 0.1, palette[subset], adj = c(0, 0))
        y <- rep(subset, each = n)
        rect(rep(0:(n-1), n), rev(y), rep(1:n, n), rev(y) - 0.5,
             col = cols[, subset], border = border)
        palette <- palette[-subset]
        cols <- cols[, -subset, drop = FALSE]
    }

    par(mfrow = c(1, 1), mar = c(5.1, 4.1, 4.1, 2.1), cex = 1)
}

hcl.swatch("qualitative")

hcl.swatch("sequential")

hcl.swatch("diverging")

hcl.swatch("divergingx")

Extracting data

We have been looking at data within the Great Barrier Reef Marine Park boundary. What if we want to extract data for some areas within that area? For example we might want to know the temperature within some of the zones within the marine park. We could crop and mask our data for each of our zones, but this could get messy if we have a lot of zones. Happily there is an easier way.

First lets load the geospatial data for two habitat protection zones in the Great Barrier Reef. This is just a small part of the Great Barrier Reef Marine Park zoning map. Surprisingly, habitat protection zones allow all types of fishing except trawling. You can view more details about the zones and the activities that are allowed here.

#load the zones data
zones <- vect("../data/gbr_habitat_protection_zones.gpkg")

#take a look
zones
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 2, 19  (geometries, attributes)
#>  extent      : 145.9373, 152.116, -22.41511, -16.38992  (xmin, xmax, ymin, ymax)
#>  source      : gbr_habitat_protection_zones.gpkg
#>  coord. ref. : lon/lat GDA2020 (EPSG:7844) 
#>  names       :        LEG_NAME          DATE SCHED_NO      SCHED_DESC PART_NO
#>  type        :           <chr>         <chr>    <int>           <chr>   <int>
#>  values      : Great Barrier ~ Gazetted 2004        1 Amalgamated Gr~       2
#>                Great Barrier ~ Gazetted 2004        1 Amalgamated Gr~       2
#>             TYPE       NAME GROUP_ID  CODE       PERMIT_DESC (and 9 more)
#>            <chr>      <chr>    <chr> <int>             <chr>             
#>  Habitat Protec~ HP-16-5126      HPZ    10   Arlington Reef~             
#>  Habitat Protec~ HP-21-5296      HPZ    10 Capricorn Channel

There are lots of columns, but lets use the “PERMIT_DESC” column, which is the zone location, to map them since it has distinct and interesting descriptions for them.

#plot the zones with our temperature data
plot(zones, "PERMIT_DESC")

Lets plot these zones over our temperature data for context.

plot(sst_cropped_masked)
lines(zones)

Coordinate reference systems

Let’s pause here for a very brief overview of coordinate reference systems!

Geographic coordinate systems: uses a 3-D surface to define locations on the Earth using longitude and latitude

Projected coordinate reference system: translates a GCS into 2-D so we can measure distances, areas, etc. in more useful units such as meters or, dare I say, feet.

https://gis.stackexchange.com/questions/149749/is-wgs84-a-coordinate-system-or-projection-system

https://r.geocompx.org/spatial-class#crs-intro

https://docs.qgis.org/3.34/en/docs/gentle_gis_introduction/coordinate_reference_systems.html

Although we have been able to plot the zones and the temperature raster together, they are actually in different coordinate reference systems (crs’s). To see what the crs of a raster or vector is we use crs().

crs(sst_cropped_masked)
#> [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

Well, that’s a lot of information. This is the “well-known text” version of the coordinate reference system. To get just some basic information that most of the time is all we need, we can do:

crs(sst_cropped_masked, describe = TRUE)
#>     name authority code area         extent
#> 1 WGS 84      EPSG 4326 <NA> NA, NA, NA, NA

crs(zones, describe = TRUE)
#>      name authority code
#> 1 GDA2020      EPSG 7844
#>                                                                                                                                                                       area
#> 1 Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore
#>                         extent
#> 1 93.41, 173.34, -8.47, -60.55

zonal(sst_cropped_masked, zones, fun = "mean")
#>   CRW_SST_2
#> 1  26.11176
#> 2  24.26957

zonal(sst_cropped_masked, zones, fun = "max")
#>   CRW_SST_2
#> 1     26.71
#> 2     24.68

zones_data_df <- extract(sst_cropped_masked, zones)
#> Warning: [extract] transforming vector data to the CRS of the raster

zones_data_df
#>     ID CRW_SST_2
#> 1    1     26.71
#> 2    1     26.61
#> 3    1     26.63
#> 4    1     26.51
#> 5    1     25.95
#> 6    1     26.16
#> 7    1     26.39
#> 8    1     26.60
#> 9    1     25.59
#> 10   1     25.79
#> 11   1     26.00
#> 12   1     26.25
#> 13   1     25.42
#> 14   1     25.63
#> 15   1     25.85
#> 16   1     26.12
#> 17   1     25.69
#> 18   2     23.97
#> 19   2     24.04
#> 20   2     24.11
#> 21   2     24.17
#> 22   2     24.21
#> 23   2     24.24
#> 24   2     24.29
#> 25   2     24.34
#> 26   2     24.39
#> 27   2     24.44
#> 28   2     23.92
#> 29   2     24.00
#> 30   2     24.06
#> 31   2     24.10
#> 32   2     24.14
#> 33   2     24.16
#> 34   2     24.21
#> 35   2     24.25
#> 36   2     24.30
#> 37   2     24.35
#> 38   2     24.40
#> 39   2     24.45
#> 40   2     23.92
#> 41   2     23.98
#> 42   2     24.02
#> 43   2     24.05
#> 44   2     24.07
#> 45   2     24.08
#> 46   2     24.12
#> 47   2     24.15
#> 48   2     24.20
#> 49   2     24.26
#> 50   2     24.32
#> 51   2     24.38
#> 52   2     24.44
#> 53   2     23.94
#> 54   2     23.97
#> 55   2     24.00
#> 56   2     24.01
#> 57   2     24.01
#> 58   2     24.01
#> 59   2     24.04
#> 60   2     24.07
#> 61   2     24.12
#> 62   2     24.18
#> 63   2     24.25
#> 64   2     24.32
#> 65   2     24.39
#> 66   2     24.45
#> 67   2     24.50
#> 68   2     23.94
#> 69   2     23.97
#> 70   2     23.98
#> 71   2     23.98
#> 72   2     23.97
#> 73   2     23.97
#> 74   2     23.99
#> 75   2     24.03
#> 76   2     24.08
#> 77   2     24.15
#> 78   2     24.22
#> 79   2     24.30
#> 80   2     24.38
#> 81   2     24.44
#> 82   2     24.48
#> 83   2     24.52
#> 84   2     24.56
#> 85   2     23.93
#> 86   2     23.97
#> 87   2     23.99
#> 88   2     23.99
#> 89   2     23.98
#> 90   2     23.99
#> 91   2     24.01
#> 92   2     24.05
#> 93   2     24.10
#> 94   2     24.18
#> 95   2     24.26
#> 96   2     24.34
#> 97   2     24.42
#> 98   2     24.47
#> 99   2     24.50
#> 100  2     24.53
#> 101  2     24.56
#> 102  2     24.58
#> 103  2     24.60
#> 104  2     23.88
#> 105  2     23.95
#> 106  2     23.99
#> 107  2     24.02
#> 108  2     24.02
#> 109  2     24.05
#> 110  2     24.08
#> 111  2     24.12
#> 112  2     24.17
#> 113  2     24.25
#> 114  2     24.34
#> 115  2     24.42
#> 116  2     24.49
#> 117  2     24.53
#> 118  2     24.55
#> 119  2     24.57
#> 120  2     24.59
#> 121  2     24.61
#> 122  2     24.62
#> 123  2     24.64
#> 124  2     23.82
#> 125  2     23.92
#> 126  2     24.00
#> 127  2     24.05
#> 128  2     24.08
#> 129  2     24.12
#> 130  2     24.16
#> 131  2     24.20
#> 132  2     24.26
#> 133  2     24.33
#> 134  2     24.42
#> 135  2     24.50
#> 136  2     24.57
#> 137  2     24.60
#> 138  2     24.61
#> 139  2     24.62
#> 140  2     24.63
#> 141  2     24.64
#> 142  2     24.65
#> 143  2     24.66
#> 144  2     24.66
#> 145  2     24.64
#> 146  2     23.77
#> 147  2     23.89
#> 148  2     24.00
#> 149  2     24.08
#> 150  2     24.13
#> 151  2     24.18
#> 152  2     24.22
#> 153  2     24.25
#> 154  2     24.31
#> 155  2     24.36
#> 156  2     24.44
#> 157  2     24.52
#> 158  2     24.58
#> 159  2     24.61
#> 160  2     24.63
#> 161  2     24.64
#> 162  2     24.65
#> 163  2     24.66
#> 164  2     24.67
#> 165  2     24.67
#> 166  2     23.72
#> 167  2     23.85
#> 168  2     23.98
#> 169  2     24.07
#> 170  2     24.13
#> 171  2     24.19
#> 172  2     24.22
#> 173  2     24.26
#> 174  2     24.30
#> 175  2     24.35
#> 176  2     24.42
#> 177  2     24.49
#> 178  2     24.54
#> 179  2     24.58
#> 180  2     24.60
#> 181  2     24.63
#> 182  2     24.65
#> 183  2     24.67
#> 184  2     24.68
#> 185  2     24.68
#> 186  2     23.67
#> 187  2     23.81
#> 188  2     23.94
#> 189  2     24.04
#> 190  2     24.10
#> 191  2     24.14
#> 192  2     24.17
#> 193  2     24.21
#> 194  2     24.25
#> 195  2     24.30
#> 196  2     24.36
#> 197  2     24.42
#> 198  2     24.47
#> 199  2     24.51
#> 200  2     24.55
#> 201  2     24.59
#> 202  2     24.63
#> 203  2     24.66
#> 204  2     24.68
#> 205  2     24.68

boxplot(CRW_SST_2 ~ ID, zones_data_df)

Rasterlayers

A very useful feature of rasters is that they can have many layers. These layers often represent different time periods, such as days or months. Lets look at a multi-layer raster which has the same temperature data that we have been looking at, but each layer represents the mean temperature in a month. We load the multi-layer raster in the same way as any other raster, using rast():

sst_monthly <- rast("../data/gbr_monthly_temp.tif")

sst_monthly
#> class       : SpatRaster 
#> dimensions  : 301, 261, 36  (nrow, ncol, nlyr)
#> resolution  : 0.05, 0.05  (x, y)
#> extent      : 142, 155.05, -25, -9.95  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : gbr_monthly_temp.tif 
#> names       : ym_202006, ym_202007, ym_202008, ym_202009, ym_202010, ym_202011, ... 
#> min values  :  19.80000,  18.51452,  19.00613,    21.096,  23.04516,    24.715, ... 
#> max values  :  29.47767,  28.65097,  28.51323,    28.864,  29.41645,    29.452, ...

We can that the raster has 36 layers. Lets plot it to see what we get.

plot(sst_monthly)

We get one map for each raster layer, but not all the data because you wouldn’t be able to see the maps if they were much smaller. We can plot a specific subset of the data using double square brackets [[]] to select a set of layers. For example, if we wanted the last 4 layers in the raster, which are layers 33 to 36:

plot(sst_monthly[[33:36]])

Happily, we can use the same functions we have learnt on this multi-layer raster. For example, we can convert all the rasters to Fahrenheit:

sst_monthly_fahrenheit <- (sst_monthly * 1.8) + 32

plot(sst_monthly_fahrenheit)

We can also crop and mask our data:

sst_monthly_cropped_masked <- crop(sst_monthly, gbr_boundary, mask = TRUE)

plot(sst_monthly_cropped_masked)

Introduction to rasters with `terra`

Jason Flower, UC Santa Barbara

2024-05-01

Prerequisites

Spatial data

Vector data

Raster data

Getting started

Exploring some real world data

Load the data

Plot the data

Load and plot some vector data

Crop and mask data

Raster values

Histogram

Frequency table

Statistics

Classifying

Raster math

Nice maps

Extracting data

Coordinate reference systems

Rasterlayers

Introduction to rasters with terra

Jason Flower, UC Santa Barbara

2024-05-01

Prerequisites

Spatial data

Vector data

Raster data

Getting started

Exploring some real world data

Load the data

Plot the data

Load and plot some vector data

Crop and mask data

Raster values

Histogram

Frequency table

Statistics

Classifying

Raster math

Nice maps

Extracting data

Coordinate reference systems

Rasterlayers

Introduction to rasters with `terra`